## Quitting Globalization: Trade-Related Job Losses, Nationalism, and Resistance
## to FDI in the United States
## Andrew Kerner, Jane Lawrence Sumner, and Yilang Feng
## Please direct all questions about R code to Jane Sumner (jlsumner@umn.edu)

library(MASS)
library(stargazer)
library(sandwich)
library(lme4)
library(aod)

####################################
###### ADDITIONAL FUNCTIONS
###### NOT BUILT INTO PACKAGES
####################################

# This function alters the transparency on the colors for some of the plots to make it
# easier to read.
# The code is from:
# http://www.dataanalytics.org.uk/Data%20Analysis/TipsAndTricks/TTR-20150531.htm
t_col <- function(color, percent = 50, name = NULL) {
  #	  color = color name
  #	percent = % transparency
  #	   name = an optional name for the color
  ## Get RGB values for named color
  rgb.val <- col2rgb(color)
  ## Make new color using input color as base and alpha set by transparency
  t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
               max = 255,
               alpha = (100-percent)*255/100,
               names = name)
  ## Save the color
  invisible(t.col)
  
}

## This function makes the rugs at the bottoms of the plots.
better.rug <- function(data,ybot,scale,lwd){
  cts <- table(data)/
    sum(table(data))
  segments(x0=as.numeric(names(cts)),
           y0=ybot,y1=ybot+cts*100/scale,lwd=lwd)
  
}
######################################################
###### DATA:
###### There are two sources of data
###### (1) Data from the Rhodium Group on 
###### Chinese investment in US, adapted from:
###### http://cim.rhg.com/
###### (used in FIgure 1)
###### (2) Data from our survey experiment
######################################################

rhodium <- read.csv("rhodium.csv",stringsAsFactors=F)

dat <- read.csv("quitting-globalization.csv",stringsAsFactors=F)

####################################
###### FIGURES AND TABLES
####################################

## Figure 1
rhodium.plot <- t(rhodium[,2:3])
colnames(rhodium.plot) <- rhodium[,1]

barplot(rhodium.plot,beside=T,ylim=c(0,50000),
        main="Chinese FDI in the United States, Greenfield and Acquisition,\n2000-2016 in USD Million",bty="n",
        col=c("gray","black"),yaxt="n")
abline(h=seq(0,50000,5000),lty=3,col="gray")
barplot(rhodium.plot,beside=T,ylim=c(0,50000),
        main="Chinese FDI in the United States, Greenfield and Acquisition,\n2000-2016 in USD Million",bty="n",
        col=c("gray","black"),yaxt="n",add=T)
legend("topleft",c("Greenfield","Acquisition"),bty="n",col=c("gray","black"),pch=15,cex=1.3,pt.cex=1.3)
axis(2,at=seq(0,50000,10000),c("$0","$10K","$20K","$30K","$40K","$50K"),las=2)


## Figure 2
hist(dat$TAApercapita,col="gray",main="Distribution of TAA Claims per Capita",freq=T,
     xlab="TAA Claims per Capita")


## Figure 3 
dat$goodforUS <- factor(dat$goodforUS,
                        levels=c("Bad","Neither","Good"))
tab <- apply(table(dat$goodforUS,dat$treatment),2,function(x){x/sum(x)})
barplot(tab,beside=T,
        main="Is This Investment Good for the US?",ylab="Proportion of Treatment Group")
legend("topleft",c("Bad","Neither","Good"),fill=gray.colors(3))

## Table 1
dat$treatment <- factor(dat$treatment,
                        levels=c("no information","foreign","Chinese"))

mod0 <- glm(investmentbad~
              treatment*TAApclogscale,data=dat,family=binomial(link="logit"))
mod1mlmcompare <- glm(investmentbad~
                        treatment*TAApclogscale+unemployedpercapita+
                        man+white+agescale+republican+manufacturing+
                        college,data=dat,family=binomial(link="logit"))

vcv0 <- vcovCL(mod0,cluster=~countyfull,data=dat)
sd0 <- sqrt(diag(vcv0))

vcv1 <- vcovCL(mod1mlmcompare,cluster=~countyfull,data=dat)
sd1 <- sqrt(diag(vcv1))

stargazer(mod0,mod1mlmcompare,se=list(sd0,sd1),no.space=T,type="html",out="table1-new.html")


## Figure 4
range.claims <- seq(min(dat$TAApclogscale,na.rm=T),max(dat$TAApclogscale,na.rm=T),.01)
vcv <- vcv1
me.foreign <- coef(mod1mlmcompare)[2]+range.claims*coef(mod1mlmcompare)[12]
varforeign <- vcv[2,2]+range.claims^2*vcv[12,12]+2*range.claims*vcv[2,12]

me.chinese <- coef(mod1mlmcompare)[3]+range.claims*coef(mod1mlmcompare)[13]
varchinese <- vcv[3,3]+range.claims^2*vcv[13,13]+2*range.claims*vcv[3,13]

plot(range.claims,me.foreign,type="l",ylim=c(-2,7),main="Marginal Effect of Treatment",
     xlab="TAA Claims per capita (logged and rescaled)",ylab="Marginal Effect of Treatment")
segments(x0=range.claims,y0=(me.foreign+1.96*sqrt(varforeign)),y1=(me.foreign-1.96*sqrt(varforeign)),
         col="gray")
lines(range.claims,me.foreign,lwd=2)
lines(range.claims,(me.foreign+1.96*sqrt(varforeign)),lty=1)
lines(range.claims,(me.foreign-1.96*sqrt(varforeign)),lty=1)
segments(x0=range.claims,y0=(me.chinese+1.96*sqrt(varchinese)),y1=(me.chinese-1.96*sqrt(varchinese)),
         col=t_col("lightgray", perc = 25, name = "lt.gray"))
lines(range.claims,me.chinese,lwd=2,lty=2)
lines(range.claims,(me.chinese+1.96*sqrt(varchinese)),lty=2)
lines(range.claims,(me.chinese-1.96*sqrt(varchinese)),lty=2)

legend("topleft",c("Foreign","Chinese"),lty=c(1,2),lwd=2)
abline(h=0,lty=3)

## Figure 5: Predicted probabilities

pp.china <- NULL
pp.foreign <- NULL
pp.domestic <- NULL
for(i in 1:1000){
  bs <- dat[sample(c(1:nrow(dat)),nrow(dat),replace=T),]
  mod.boot <- glm(investmentbad~
                    treatment*TAApclogscale+unemployedpercapita+
                    man+white+agescale+republican+manufacturing+
                    college,data=bs,family=binomial(link="logit"))

  range.claims <- seq(min(dat$TAApclogscale,na.rm=T),max(dat$TAApclogscale,na.rm=T),.01)
  
  hyp.domestic <- cbind(1,
                        0,0,
                        range.claims,
                        mean(dat$unemployedpercapita,na.rm=T),
                        mean(dat$man,na.rm=T),
                        mean(dat$white,na.rm=T),
                        mean(dat$agescale,na.rm=T),
                        mean(dat$republican,na.rm=T),
                        mean(dat$manufacturing,na.rm=T),
                        mean(dat$college,na.rm=T),
                        0,0)
  xb.domestic <- hyp.domestic %*% coef(mod.boot)
  hyp.chinese <- cbind(1,
                       0,1,
                       range.claims,
                       mean(dat$unemployedpercapita,na.rm=T),
                       mean(dat$man,na.rm=T),
                       mean(dat$white,na.rm=T),
                       mean(dat$agescale,na.rm=T),
                       mean(dat$republican,na.rm=T),
                       mean(dat$manufacturing,na.rm=T),
                       mean(dat$college,na.rm=T),
                       0,range.claims)
  xb.chinese <- hyp.chinese %*% coef(mod.boot)
  hyp.foreign <- cbind(1,
                       1,0,
                       range.claims,
                       mean(dat$unemployedpercapita,na.rm=T),
                       mean(dat$man,na.rm=T),
                       mean(dat$white,na.rm=T),
                       mean(dat$agescale,na.rm=T),
                       mean(dat$republican,na.rm=T),
                       mean(dat$manufacturing,na.rm=T),
                       mean(dat$college,na.rm=T),
                       range.claims,0)
  xb.foreign <- hyp.foreign %*% coef(mod.boot)
  
  pp.china <- cbind(pp.china,plogis(xb.chinese))
  pp.foreign <- cbind(pp.foreign,plogis(xb.foreign))
  pp.domestic <- cbind(pp.domestic,plogis(xb.domestic))
}

hyp.domestic <- cbind(1,
                      0,0,
                      range.claims,
                      mean(dat$unemployedpercapita,na.rm=T),
                      mean(dat$man,na.rm=T),
                      mean(dat$white,na.rm=T),
                      mean(dat$agescale,na.rm=T),
                      mean(dat$republican,na.rm=T),
                      mean(dat$manufacturing,na.rm=T),
                      mean(dat$college,na.rm=T),
                      0,0)
xb.domestic <- hyp.domestic %*% coef(mod1mlmcompare)
hyp.chinese <- cbind(1,
                     0,1,
                     range.claims,
                     mean(dat$unemployedpercapita,na.rm=T),
                     mean(dat$man,na.rm=T),
                     mean(dat$white,na.rm=T),
                     mean(dat$agescale,na.rm=T),
                     mean(dat$republican,na.rm=T),
                     mean(dat$manufacturing,na.rm=T),
                     mean(dat$college,na.rm=T),
                     0,range.claims)
xb.chinese <- hyp.chinese %*% coef(mod1mlmcompare)
hyp.foreign <- cbind(1,
                     1,0,
                     range.claims,
                     mean(dat$unemployedpercapita,na.rm=T),
                     mean(dat$man,na.rm=T),
                     mean(dat$white,na.rm=T),
                     mean(dat$agescale,na.rm=T),
                     mean(dat$republican,na.rm=T),
                     mean(dat$manufacturing,na.rm=T),
                     mean(dat$college,na.rm=T),
                     range.claims,0)
xb.foreign <- hyp.foreign %*% coef(mod1mlmcompare)


lowerCI.domestic <- apply(pp.domestic,1,quantile,.025,na.rm=T)
upperCI.domestic <- apply(pp.domestic,1,quantile,.975,na.rm=T)
lowerCI.china <- apply(pp.china,1,quantile,.025,na.rm=T)
upperCI.china <- apply(pp.china,1,quantile,.975,na.rm=T)
lowerCI.foreign <- apply(pp.foreign,1,quantile,.025,na.rm=T)
upperCI.foreign <- apply(pp.foreign,1,quantile,.975,na.rm=T)

plot(range.claims,plogis(xb.chinese),ylim=c(0,.8),type="l",lwd=2,
     xlab="range of logged TAA claims",ylab="Probability of Believing Investment is Bad",
     main="Predicted Probabilities\nby Treatment Group")
segments(x0=range.claims[seq(1,length(range.claims),10)],
         y0=lowerCI.china[seq(1,length(range.claims),10)],
         y1=upperCI.china[seq(1,length(range.claims),10)])
lines(range.claims,plogis(xb.domestic),lwd=2,lty=2)
segments(x0=range.claims[seq(1,length(range.claims),10)],
         y0=lowerCI.domestic[seq(1,length(range.claims),10)],
         y1=upperCI.domestic[seq(1,length(range.claims),10)],lty=2)
lines(range.claims,plogis(xb.foreign),lwd=3,lty=2)
segments(x0=range.claims[seq(1,length(range.claims),10)],
         y0=lowerCI.foreign[seq(1,length(range.claims),10)],
         y1=upperCI.foreign[seq(1,length(range.claims),10)],lty=3)

legend("topright",lwd=2,lty=c(1,2,3),c("China Treatment","Foreign Treatment","No Info Treatment"))

## Table 2
mod.culture.nat <- glm(investmentbad~
                         treatment*TAApclogscale+unemployedpercapita+
                         man+white+age+republican+manufacturing+
                         college,data=dat[which(dat$culturalnationalism==1),],family=binomial(link="logit"))

vcv1 <- vcovCL(mod.culture.nat,cluster=~countyfull,data=dat)
sd1 <- sqrt(diag(vcv1))

mod.culture.nat2 <- glm(investmentbad~
                          treatment*TAApclogscale,data=dat[which(dat$culturalnationalism==1),],family=binomial(link="logit"))

mod.culture.nn <- glm(investmentbad~
                        treatment*TAApclogscale+unemployedpercapita+
                        man+white+age+republican+manufacturing+
                        college,data=dat[which(dat$culturalnationalism==0),],family=binomial(link="logit"))
mod.culture.nn2 <- glm(investmentbad~
                         treatment*TAApclogscale,data=dat[which(dat$culturalnationalism==0),],family=binomial(link="logit"))
vcv2 <- vcovCL(mod.culture.nn,cluster=~countyfull,data=dat)
sd2 <- sqrt(diag(vcv2))

stargazer(mod.culture.nat,mod.culture.nn,
          se=list(sd1,sd2),type="html",no.space=T,out="apsa-table-culture2.html")

## Figure 6
## (In hindsight, I should have written a function to do this. But here we are.)
vcv <- vcv1
mod0 <- mod.culture.nat
me.foreign <- coef(mod0)[2]+range.claims*coef(mod0)[12]
varforeign <- vcv[2,2]+range.claims^2*vcv[12,12]+2*range.claims*vcv[2,12]

me.chinese <- coef(mod0)[3]+range.claims*coef(mod0)[13]
varchinese <- vcv[3,3]+range.claims^2*vcv[13,13]+2*range.claims*vcv[3,13]

par(mfrow=c(2,2))
plot(range.claims,me.foreign,type="l",
     main="Foreign Treatment\nNationalists",
     xlab="TAA Claims per capita (logged and rescaled)",ylab="Marginal Effect of Treatment",
     yaxs="i",bty="n",ylim=c(-5,9))
segments(x0=range.claims,y0=(me.foreign+1.96*sqrt(varforeign)),y1=(me.foreign-1.96*sqrt(varforeign)),
         col="gray")
lines(range.claims,me.foreign,lwd=2)
lines(range.claims,(me.foreign+1.96*sqrt(varforeign)),lty=1)
lines(range.claims,(me.foreign-1.96*sqrt(varforeign)),lty=1)
abline(h=0,lty=3)
better.rug(data=dat$TAApclogscale[which(dat$treatment=="foreign" &
                                          dat$culturalnationalism==1)],
           ybot=-5,scale=4,lwd=4)

plot(range.claims,me.chinese,type="l",ylim=c(-5,9),
     main="Chinese Treatment\nNationalists",
     xlab="TAA Claims per capita (logged and rescaled)",ylab="Marginal Effect of Treatment",
     yaxs="i",bty="n")
segments(x0=range.claims,y0=(me.chinese+1.96*sqrt(varchinese)),y1=(me.chinese-1.96*sqrt(varchinese)),
         col="gray")
lines(range.claims,me.chinese,lwd=2,lty=1)
lines(range.claims,(me.chinese+1.96*sqrt(varchinese)),lty=1)
lines(range.claims,(me.chinese-1.96*sqrt(varchinese)),lty=1)
abline(h=0,lty=3)
better.rug(data=dat$TAApclogscale[which(dat$treatment=="Chinese" &
                                          dat$culturalnationalism==1)],
           ybot=-5,scale=4,lwd=4)

vcv <- vcv2
mod0 <- mod.culture.nn
me.foreign <- coef(mod0)[2]+range.claims*coef(mod0)[12]
varforeign <- vcv[2,2]+range.claims^2*vcv[12,12]+2*range.claims*vcv[2,12]

me.chinese <- coef(mod0)[3]+range.claims*coef(mod0)[13]
varchinese <- vcv[3,3]+range.claims^2*vcv[13,13]+2*range.claims*vcv[3,13]

plot(range.claims,me.foreign,type="l",ylim=c(-5,9),main="Foreign Treatment\nNon-Nationalists",
     xlab="TAA Claims per capita (logged and rescaled)",ylab="Marginal Effect of Treatment",
     xaxs="i",bty="n",yaxs="i")
segments(x0=range.claims,y0=(me.chinese+1.96*sqrt(varchinese)),
         y1=(me.chinese-1.96*sqrt(varchinese)),
         col="gray")
lines(range.claims,me.chinese,lwd=2,lty=1)
lines(range.claims,(me.chinese+1.96*sqrt(varchinese)),lty=1)
lines(range.claims,(me.chinese-1.96*sqrt(varchinese)),lty=1)
abline(h=0,lty=3)
better.rug(data=dat$TAApclogscale[which(dat$treatment=="foreign" &
                                          dat$culturalnationalism==0)],
           ybot=-5,scale=4,lwd=4)

plot(range.claims,me.chinese,type="l",ylim=c(-5,9),main="Chinese Treatment\nNon-Nationalists",
     xlab="TAA Claims per capita (logged and rescaled)",ylab="Marginal Effect of Treatment",
     xaxs="i",bty="n",yaxs="i")
segments(x0=range.claims,y0=(me.foreign+1.96*sqrt(varforeign)),y1=(me.foreign-1.96*sqrt(varforeign)),
         col="gray")
lines(range.claims,me.foreign,lwd=2)
lines(range.claims,(me.foreign+1.96*sqrt(varforeign)),lty=1)
lines(range.claims,(me.foreign-1.96*sqrt(varforeign)),lty=1)
abline(h=0,lty=3)
better.rug(data=dat$TAApclogscale[which(dat$treatment=="Chinese" &
                                          dat$culturalnationalism==0)],
           ybot=-5,scale=4,lwd=4)
par(mfrow=c(1,1))

## Figure 7: predicted probabilities
## Note: To make the legend not overlap the CI in the top right panel, I use windows()
## because I am a PC user and it allows me to stretch the window to the right size.
## I understand this isn't exactly the same for Macs, so your mileage may vary.

pp.china <- NULL
pp.foreign <- NULL
pp.domestic <- NULL
for(i in 1:1000){
  # nationalists
  subdat <- dat[which(dat$culturalnationalism==1),]
  bs <- subdat[sample(c(1:nrow(subdat)),nrow(subdat),replace=T),]
  mod.boot <- glm(investmentbad~
                    treatment*TAApclogscale+unemployedpercapita+
                    man+white+age+republican+manufacturing+
                    college,
                  data=bs,family=binomial(link="logit"))
  
  
  hyp.domestic <- cbind(1,
                        0,0,
                        range.claims,
                        mean(dat$unemployedpercapita,na.rm=T),
                        mean(dat$man,na.rm=T),
                        mean(dat$white,na.rm=T),
                        mean(dat$age,na.rm=T),
                        mean(dat$republican,na.rm=T),
                        mean(dat$manufacturing,na.rm=T),
                        mean(dat$college,na.rm=T),
                        0,0)
  xb.domestic <- hyp.domestic %*% coef(mod.boot)
  hyp.chinese <- cbind(1,
                       0,1,
                       range.claims,
                       mean(dat$unemployedpercapita,na.rm=T),
                       mean(dat$man,na.rm=T),
                       mean(dat$white,na.rm=T),
                       mean(dat$age,na.rm=T),
                       mean(dat$republican,na.rm=T),
                       mean(dat$manufacturing,na.rm=T),
                       mean(dat$college,na.rm=T),
                       0,range.claims)
  xb.chinese <- hyp.chinese %*% coef(mod.boot)
  hyp.foreign <- cbind(1,
                       1,0,
                       range.claims,
                       mean(dat$unemployedpercapita,na.rm=T),
                       mean(dat$man,na.rm=T),
                       mean(dat$white,na.rm=T),
                       mean(dat$age,na.rm=T),
                       mean(dat$republican,na.rm=T),
                       mean(dat$manufacturing,na.rm=T),
                       mean(dat$college,na.rm=T),
                       range.claims,0)
  xb.foreign <- hyp.foreign %*% coef(mod.boot)
  
  
  pp.china <- cbind(pp.china,plogis(xb.chinese))
  pp.foreign <- cbind(pp.foreign,plogis(xb.foreign))
  pp.domestic <- cbind(pp.domestic,plogis(xb.domestic))
}

hyp.domestic <- cbind(1,
                      0,0,
                      range.claims,
                      mean(dat$unemployedpercapita,na.rm=T),
                      mean(dat$man,na.rm=T),
                      mean(dat$white,na.rm=T),
                      mean(dat$age,na.rm=T),
                      mean(dat$republican,na.rm=T),
                      mean(dat$manufacturing,na.rm=T),
                      mean(dat$college,na.rm=T),
                      0,0)
xb.domestic <- hyp.domestic %*% coef(mod.culture.nat)
hyp.chinese <- cbind(1,
                     0,1,
                     range.claims,
                     mean(dat$unemployedpercapita,na.rm=T),
                     mean(dat$man,na.rm=T),
                     mean(dat$white,na.rm=T),
                     mean(dat$age,na.rm=T),
                     mean(dat$republican,na.rm=T),
                     mean(dat$manufacturing,na.rm=T),
                     mean(dat$college,na.rm=T),
                     0,range.claims)
xb.chinese <- hyp.chinese %*% coef(mod.culture.nat)
hyp.foreign <- cbind(1,
                     1,0,
                     range.claims,
                     mean(dat$unemployedpercapita,na.rm=T),
                     mean(dat$man,na.rm=T),
                     mean(dat$white,na.rm=T),
                     mean(dat$age,na.rm=T),
                     mean(dat$republican,na.rm=T),
                     mean(dat$manufacturing,na.rm=T),
                     mean(dat$college,na.rm=T),
                     range.claims,0)
xb.foreign <- hyp.foreign %*% coef(mod.culture.nat)


lowerCI.domestic <- apply(pp.domestic,1,quantile,.025)
upperCI.domestic <- apply(pp.domestic,1,quantile,.975)
lowerCI.china <- apply(pp.china,1,quantile,.025)
upperCI.china <- apply(pp.china,1,quantile,.975)
lowerCI.foreign <- apply(pp.foreign,1,quantile,.025)
upperCI.foreign <- apply(pp.foreign,1,quantile,.975)


pp.china2 <- NULL
pp.foreign2 <- NULL
pp.domestic2 <- NULL
for(i in 1:1000){
  subdat <- dat[which(dat$culturalnationalism==0),]
  bs <- subdat[sample(c(1:nrow(subdat)),nrow(subdat),replace=T),]
  mod.boot <- glm(investmentbad~
                    treatment*TAApclogscale+unemployedpercapita+
                    man+white+age+republican+manufacturing+
                    college,
                  data=bs,family=binomial(link="logit"))
  
  
  hyp.domestic <- cbind(1,
                        0,0,
                        range.claims,
                        mean(dat$unemployedpercapita,na.rm=T),
                        mean(dat$man,na.rm=T),
                        mean(dat$white,na.rm=T),
                        mean(dat$age,na.rm=T),
                        mean(dat$republican,na.rm=T),
                        mean(dat$manufacturing,na.rm=T),
                        mean(dat$college,na.rm=T),
                        0,0)
  xb.domestic2 <- hyp.domestic %*% coef(mod.boot)
  hyp.chinese <- cbind(1,
                       0,1,
                       range.claims,
                       mean(dat$unemployedpercapita,na.rm=T),
                       mean(dat$man,na.rm=T),
                       mean(dat$white,na.rm=T),
                       mean(dat$age,na.rm=T),
                       mean(dat$republican,na.rm=T),
                       mean(dat$manufacturing,na.rm=T),
                       mean(dat$college,na.rm=T),
                       0,range.claims)
  xb.chinese2 <- hyp.chinese %*% coef(mod.boot)
  hyp.foreign <- cbind(1,
                       1,0,
                       range.claims,
                       mean(dat$unemployedpercapita,na.rm=T),
                       mean(dat$man,na.rm=T),
                       mean(dat$white,na.rm=T),
                       mean(dat$age,na.rm=T),
                       mean(dat$republican,na.rm=T),
                       mean(dat$manufacturing,na.rm=T),
                       mean(dat$college,na.rm=T),
                       range.claims,0)
  xb.foreign2 <- hyp.foreign %*% coef(mod.boot)
  
  pp.china2 <- cbind(pp.china2,plogis(xb.chinese2))
  pp.foreign2 <- cbind(pp.foreign2,plogis(xb.foreign2))
  pp.domestic2 <- cbind(pp.domestic2,plogis(xb.domestic2))
}

hyp.domestic <- cbind(1,
                      0,0,
                      range.claims,
                      mean(dat$unemployedpercapita,na.rm=T),
                      mean(dat$man,na.rm=T),
                      mean(dat$white,na.rm=T),
                      mean(dat$age,na.rm=T),
                      mean(dat$republican,na.rm=T),
                      mean(dat$manufacturing,na.rm=T),
                      mean(dat$college,na.rm=T),
                      0,0)
xb.domestic2 <- hyp.domestic %*% coef(mod.culture.nn)
hyp.chinese <- cbind(1,
                     0,1,range.claims,
                     mean(dat$unemployedpercapita,na.rm=T),
                     mean(dat$man,na.rm=T),
                     mean(dat$white,na.rm=T),
                     mean(dat$age,na.rm=T),
                     mean(dat$republican,na.rm=T),
                     mean(dat$manufacturing,na.rm=T),
                     mean(dat$college,na.rm=T),
                     0,range.claims)
xb.chinese2 <- hyp.chinese %*% coef(mod.culture.nn)
hyp.foreign <- cbind(1,
                     1,0,range.claims,
                     mean(dat$unemployedpercapita,na.rm=T),
                     mean(dat$man,na.rm=T),
                     mean(dat$white,na.rm=T),
                     mean(dat$age,na.rm=T),
                     mean(dat$republican,na.rm=T),
                     mean(dat$manufacturing,na.rm=T),
                     mean(dat$college,na.rm=T),
                     range.claims,0)
xb.foreign2 <- hyp.foreign %*% coef(mod.culture.nn)


lowerCI.domestic2 <- apply(pp.domestic2,1,quantile,.025)
upperCI.domestic2 <- apply(pp.domestic2,1,quantile,.975)
lowerCI.china2 <- apply(pp.china2,1,quantile,.025)
upperCI.china2 <- apply(pp.china2,1,quantile,.975)
lowerCI.foreign2 <- apply(pp.foreign2,1,quantile,.025)
upperCI.foreign2 <- apply(pp.foreign2,1,quantile,.975)

# 2 denotes non-nationalists
par(mfrow=c(2,2))
plot(range.claims,plogis(xb.foreign),ylim=c(0,1),type="l",lwd=2,
     xlab="TAA Claims per capita (logged and rescaled)"
     ,ylab="Probability of Believing Investment is Bad",
     main="Predicted Probabilities for Nationalists\nForeign v. Generic")
segments(x0=range.claims,
         y0=lowerCI.foreign,
         y1=upperCI.foreign,col="lightgray")
lines(range.claims,plogis(xb.foreign),lwd=2)
lines(range.claims,lowerCI.foreign,lty=1)
lines(range.claims,upperCI.foreign,lty=1)
segments(x0=range.claims,
         y0=lowerCI.domestic,
         y1=upperCI.domestic,col="gray")
lines(range.claims,plogis(xb.domestic),lwd=2,lty=2)
lines(range.claims,lowerCI.domestic,lty=1)
lines(range.claims,upperCI.domestic,lty=1)
legend("topright",lty=c(1,2),lwd=2,c("Foreign","Generic"))

plot(range.claims,plogis(xb.chinese),ylim=c(0,1),type="l",lwd=2,
     xlab="TAA Claims per capita (logged and rescaled)"
     ,ylab="Probability of Believing Investment is Bad",
     main="Predicted Probabilities for Nationalists\n Chinese v. Generic")
segments(x0=range.claims,
         y0=lowerCI.china,
         y1=upperCI.china,col="lightgray")
lines(range.claims,plogis(xb.chinese),lwd=2)
lines(range.claims,lowerCI.china,lty=1)
lines(range.claims,upperCI.china,lty=1)
segments(x0=range.claims,
         y0=lowerCI.domestic,
         y1=upperCI.domestic,col="gray")
lines(range.claims,plogis(xb.domestic),lwd=2,lty=2)
lines(range.claims,lowerCI.domestic,lty=1)
lines(range.claims,upperCI.domestic,lty=1)
legend("topright",lty=c(1,2),lwd=2,c("Chinese","Generic"))

plot(range.claims,plogis(xb.foreign2),ylim=c(0,1),type="l",lwd=2,
     xlab="TAA Claims per capita (logged and rescaled)"
     ,ylab="Probability of Believing Investment is Bad",
     main="Predicted Probabilities for Non-Nationalists\nForeign v. Generic",xasx="i")
segments(x0=range.claims,
         y0=lowerCI.foreign2,
         y1=upperCI.foreign2,col="lightgray")
lines(range.claims,plogis(xb.foreign2),lwd=2)
lines(range.claims,lowerCI.foreign2,lty=1)
lines(range.claims,upperCI.foreign2,lty=1)
segments(x0=range.claims,
         y0=lowerCI.domestic2,
         y1=upperCI.domestic2,col="gray")
lines(range.claims,plogis(xb.domestic2),lwd=2,lty=2)
lines(range.claims,lowerCI.domestic2,lty=1)
lines(range.claims,upperCI.domestic2,lty=1)
legend("topright",lty=c(1,2),lwd=2,c("Foreign","Generic"))

plot(range.claims,plogis(xb.chinese2),ylim=c(0,1),type="l",lwd=2,
     xlab="TAA Claims per capita (logged and rescaled)"
     ,ylab="Probability of Believing Investment is Bad",
     main="Predicted Probabilities for Non-Nationalists\n Chinese v. Generic")
segments(x0=range.claims,
         y0=lowerCI.china2,
         y1=upperCI.china2,col="lightgray")
lines(range.claims,plogis(xb.chinese2),lwd=2)
lines(range.claims,lowerCI.china2,lty=1)
lines(range.claims,upperCI.china2,lty=1)
segments(x0=range.claims,
         y0=lowerCI.domestic2,
         y1=upperCI.domestic2,col="gray")
lines(range.claims,plogis(xb.domestic2),lwd=2,lty=2)
lines(range.claims,lowerCI.domestic2,lty=1)
lines(range.claims,upperCI.domestic2,lty=1)
legend("topright",lty=c(1,2),lwd=2,c("Chinese","Generic"))
par(mfrow=c(1,1))

####### Appendix A: Balance
## Table A1
bal1 <- glm(dat$treatment=="foreign"~TAApclogscale+unemployedpercapita+
              man+white+age+republican+manufacturing+
              college,
            family=binomial(link="logit"),data=dat)
bal2 <- glm(dat$treatment=="Chinese"~TAApclogscale+unemployedpercapita+
              man+white+age+republican+manufacturing+
              college,
            family=binomial(link="logit"),data=dat)
bal3 <- glm(dat$treatment=="no information"~TAApclogscale+unemployedpercapita+
              man+white+age+republican+manufacturing+
              college,
            family=binomial(link="logit"),data=dat)

stargazer(bal1,bal2,bal3,no.space=T,type="html",out="balance.html")

## Table A2
variables <- cbind(dat$TAApclogscale,dat$unemployedpercapita,
                   dat$man,dat$white,dat$agescale,dat$republican,dat$manufacturing,
                   dat$college)
p <- matrix(NA,ncol=3,nrow=8)
for(i in 1:ncol(variables)){
  p[i,] <- c(t.test(variables[,i][which(dat$treatment=="foreign")],
                    variables[,i][which(dat$treatment=="Chinese")])$p.value,
             t.test(variables[,i][which(dat$treatment=="foreign")],
                    variables[,i][which(dat$treatment=="no information")])$p.value,
             t.test(variables[,i][which(dat$treatment=="no information")],
                    variables[,i][which(dat$treatment=="Chinese")])$p.value)
  
}
p

## Figure A1
plot(density(dat$TAApclogscale[which(dat$treatment=="foreign")],na.rm=T),ylim=c(0,.6),
     main="TAA Claims per capita (logged, rescaled)",xlab="TAA Claims per capita (logged, rescaled)")
abline(v=mean(dat$TAApclogscale[which(dat$treatment=="foreign")],na.rm=T))
lines(density(dat$TAApclogscale[which(dat$treatment=="no information")],na.rm=T),lty=3,lwd=2)
abline(v=mean(dat$TAApclogscale[which(dat$treatment=="no information")],na.rm=T),lty=3,lwd=2)
legend("topright",lty=c(1,3),lwd=c(1,2),c("Foreign","No Information"),bty="n")


####### Appendix B: Robustness Checks
## Table B1: model from body of paper with SE clustering at different levels
vcv.unclustered <- vcov(mod1mlmcompare)
sd.unclustered <- sqrt(diag(vcv.unclustered))
vcv.county <- vcovCL(mod1mlmcompare,cluster=~countyfull,data=dat)
sd.county <- sqrt(diag(vcv.county))
vcv.state <- vcovCL(mod1mlmcompare,cluster=~state,data=dat)
sd.state <- sqrt(diag(vcv.state))
vcv.treatment <- vcovCL(mod1mlmcompare,cluster=~treatment,data=dat)
sd.treatment <- sqrt(diag(vcv.treatment))

stargazer(mod1mlmcompare,mod1mlmcompare,mod1mlmcompare,mod1mlmcompare,
          se=list(sd.unclustered,
                  sd.county,
                  sd.state,
                  sd.treatment),no.space=T,type="html",out="clustering.html")

## Table B2: baseline model  v random intercepts v fixed effects/dummies
# This is the model from the body of the paper.
mod1mlmcompare <- glm(investmentbad~
                        treatment*TAApclogscale+unemployedpercapita+
                        man+white+agescale+republican+manufacturing+
                        college,data=dat,family=binomial(link="logit"))

mod1mlm <- glmer(investmentbad~
                   treatment*TAApclogscale+unemployedpercapita+
                   man+white+agescale+republican+manufacturing+
                   college+(1|state),data=dat,family=binomial(link="logit"))

mod1stateFE <- glm(investmentbad~
                     treatment*TAApclogscale+unemployedpercapita+
                     man+white+agescale+republican+manufacturing+
                     college+as.factor(state),data=dat,family=binomial(link="logit"))

# The model below is commented out as evidence that it does not converge.
#mod1countyFE <- glm(investmentbad~
#                      treatment*TAApclogscale+unemployedpercapita+
#                      man+white+agescale+republican+manufacturing+
#                      college+as.factor(countyfull),data=dat,family=binomial(link="logit"))

mod1countyFE2 <- glm(investmentbad~
                       treatment*TAApclogscale+unemployedpercapita+
                       man+white+agescale+republican+manufacturing+
                       college+
                       miami+lac+cook+sdc+nyc+maricopa+clark,data=dat,
                       family=binomial(link="logit"))


stargazer(mod1mlmcompare,mod1mlm,mod1stateFE,mod1countyFE2,no.space=T,type="html",
          out="robustness1.html")

####### Appendix C: three-way interactions
## Table C1
mod.culture.int <- glm(investmentbad~
                         treatment*TAApclogscale*culturalnationalism+unemployedpercapita+
                         man+white+age+republican+manufacturing+
                         college,data=dat,family=binomial(link="logit"))
vcv1 <- vcovCL(mod.culture.int,cluster=~countyfull,data=dat)
sd1 <- sqrt(diag(vcv1))

wald.test(Sigma=vcv1,b=coef(mod.culture.int),Terms=c(19,14))
wald.test(Sigma=vcv1,b=coef(mod.culture.int),Terms=c(18,15))
wald.test(Sigma=vcv1,b=coef(mod.culture.int),Terms=c(17,4))

stargazer(mod.culture.int,se=list(sd1),type="html",no.space=T,out="apsa-table-3way.html")


## end of file.